Exact treatment of exciton-polaron formation by Diagrammatic Monte Carlo 
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We develop an approximation-free Diagrammatic Monte Carlo technique to study fermionic par- 
ticles interacting with each other simultaneously through both an attractive Coulomb potential 
and bosonic excitations of the underlying medium. Exemplarily we apply the method to the 
long-standing exciton-polaron problem and present numerically exact results for the wave func- 
tion, ground-state energy, binding energy and effective mass of this quasiparticle. Focusing on the 
electron-hole pair bound-state formation, we discuss various limiting cases of a generic exciton- 
polaron model. The frequently used instantaneous approximation to the retarded interaction due to 
the phonon exchange is found to be of very limited applicability. For the case of a light electron and 
heavy hole the system is well approximated by a particle in the field of a static attractive impurity. 
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PACS numbers: 71.35.-y, 02.70.Ss, 71.38.-k 

The problem of two quasiparticles (QPs) interacting 
via instantaneous Coulomb interaction and a retarded 
exchange of bosons is a tremendously difficult stum- 
bling block in solid state many-particle physics [2, 0, Q • 
The variety of such objects range from exciton-polarons 
(EX-P) in semiconductors, where (opposite) charged 
holes (H) and electrons (E) couple to lattice vibrations 
(phonons) 0, to more exotic situations e.g. in the con- 
text of high-Tc superconductivity, where two holes may 
form a bound state in an antiferromagnetically corre- 
lated background due to exchange of spin excitations 
(magnons) [4|. 

First attempts to tackle the EX-P problem are re- 
stricted to low-dimensional cases, and reduce the two 
QPs to a preformed structureless QP object 0, 
Moreover the rather crude adiabatic approximation was 
frequently used 0], as well as simple variational ap- 
proaches [fjj]. Quite recently a quantum Monte Carlo 
study has been carried out [7| , but for a one-dimensional 
(ID) model with simplified Coulomb and E/H- phonon 
interactions. In any case, a general approximation- 
free method for treating a system of interacting QPs in 
bosonic fields is missing. Even the interaction of two 
QPs through an instantaneous potential creates enor- 
mous technical difficulties. For 3D situations with realis- 
tic QP dispersions, at the moment, Diagrammatic Monte 
Carlo (DMC) [H and Bethe-Salpeter methods seem to 
be the most promising techniques to address such prob- 
lems. Within a Bethe-Salpeter based approach coupling 
to phonons can be introduced only in a less controllable 
phenomenological way [10j and we are not aware of any 
generalization to a true full-scale problem. By contrast, 
we will demonstrate that a corresponding generalization 
of the DMC method can be done in a rigorous way. 

In the present paper we develop such a general 
approximation-free DMC technique and apply it for the 
first time to the highly non-trivial EX-P problem. Note 



that different from direct-space DMC our method 
is realized in momentum space and will thus not be re- 
stricted to the treatment of finite systems. Furthermore 
the proposed momentum space DMC approach is capable 
of describing dispersive fermions and bosons, as well as 
(long-range) interactions between those objects in any di- 
mension. It is used to obtain the EX-P wave function, en- 
ergy, and mass, also for the case when an E-H bound state 
arises due to cooperative effect of short-range Coulomb 
attraction and exchange of phonons. In addition we de- 
termine the ground-state phase diagram for a restricted 
EX-P model with contact Coulomb and particle-phonon 
interaction, discuss important limiting cases and the va- 
lidity of approximative solutions. 

To this end let us start from the following Hamiltonian 
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Here ek (hk) annihilates an E (H) in the conduction (va- 
lence) band with dispersion e c (k) (£„(k)) and & q is the 
corresponding annihilation operator for a phonon with 
momentum q. In (JT]) , U (p, k, k') describes the attractive 
interband interaction and <? e (q) (ffh(q)) the E(H)-phonon 
coupling. N denotes the number of lattice sites. We work 
in the thermodynamic limit N — > oo. 

An energy (^(p)) momentum (p) eigenstate (\v;p)) 
of H can be expressed as linear combination 
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of E-H pair basis states in the presence of M 
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phonons (having momenta qi, . . . , qi/), where Y„ 
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1 qj llj=i °qj- 



M{q} 
pk 



The M = term is under- 

ek^-p-k- Within this basis we de- 
fine the two-particle imaginary-time Green function (GF) 
with center-of-mass momentum p as 
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where Y(t) — e Ye~ (r > 0), and |vac) is a di- 
rect product of the phonon vacuum and completely filled 
(empty) valence (conduction) bands. Rewriting (J3|) in 
interaction representation and expanding it in terms of 
both Coulomb interaction U and E(H)-phonon couplings 
g e .hi one arrives at a series of phonon-dressed ladder-type 
Feynman diagrams cf. Fig.[T] The weight attributed 
to a given diagram is the product of the interaction ver- 
tices (£/(p,k, k'), ge(q), and g?i(q)) and Matsubara GFs 
of holes, electrons, and phonons with the corresponding 
imaginary times and momenta subjected to momentum 
conservation imposed by the Hamiltonian |T]). 




FIG. 1: A typical diagram for G^°(r). Solid (dashed) line 
represent E (H) propagators, solid circles (squares) designate 
Coulomb (QP-phonon) interactions, and dotted lines are the 
phonon propagators. Imaginary time runs from left to right. 



momentum of E, H, and phonons is equal to the center- 
of-mass momentum p. 






FIG. 2: The way how momentum conservation is ensured (b- 
c) when a nondiagonal phonon is added to configuration (a). 
Diagrams with a non-diagonal phonon line and no Coulomb 
vertices are prohibited by momentum conservation in 

Now the EX-P's energy, effective mass, and wave func- 
tion can be found by DMC sampling of GFs at times 
larger than the reciprocal energy difference between the 
ground |g.s.) and the first excited state, r > Ti; m . Insert- 
ing the complete set ([2]) into ([3]), we have 
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In the numerical work, the Monte Carlo updates of 
diagonal phonon propagators (connecting two points T\ 
and T2 of the same QP propagator) were performed by the 
DMC technique developed for a polaron 

[H, El, while 

Coulomb vertices are updated as for the pure exciton 
problem 0. The new update for nondiagonal phonon 
propagators (see Fig. [5]) requires special care to maintain 
the momentum conservation. When updating diagonal 
phonon contributions, momentum conservation is simply 
achieved by subtracting the phonon momentum from all 
QP propagators between n and T2 3, 14 1. Such strat- 
egy is not suitable for the non-diagonal phonon propa- 
gators since the phonon momentum is taken from one 
QP line and absorbed by another one. The problem can 
be solved, however, by absorbing the momentum transfer 
into the Coulomb vertex which, in circular representation 
for the GF [H, [l4| , always appears either to the left of the 
phonon propagator (r a in Fig. [5^) or between t\ and T2 
(re, in Fig. [2^i). As illustrated in Fig.[2j3, in the first case 
the phonon momentum Q is subtracted (added) to (from) 
the E (H) propagators located between the Coulomb ver- 
tex at r a and n (T2). In the second case, the E (H) 
momenta are changed as k — > k — Q for r' G [t^ , T2] 
(t' € [ti,t;,]), see Fig. Note that at any t' the total 



For r > Tii m the GF projects onto the ground state in 
the p-sector [H [l3j : 
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Due to the normalization X)m=o S{ q }k £ 
the sum of all possible M-phonon GFs [l 
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has an especially simple asymptotic form, <3 p (t — > 00) — > 
e --rSg. s .(p) According to ([5]) and ((6|) the estimators for 
the amplitudes £pk^(g- s -) are related to the distribu- 
tion of variables M, {q} and k which are all gener- 
ated by the DMC algorithm Gp^ q} (r)/0 p (r) \ T ^oo^> 

2 

£pk {q} fe- s -) • Since the whole set {^pk {q} (g- s -)}, defin- 
ing the inner structure of the EX-P state, is difficult to 

visualize, we introduce integrated quantities. For exam- 
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pie, the integrated Z-factor Z p M ^ = Sk{ q } Cpk {q} (g- S -) 
measures the partial weights of M-phonon configurations 
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FIG. 3: (Color online) Binding energy £ex-p (in units of the 
gap E g ) and effective mass mEx-p (in units of mo = E CtV /6) 
as functions of E/H-phonon coupling A for SI = 0.25E g . The 
inset shows the EX-P energy -E g . s . (circles) compared to 2Ehp 
(squares) at U/E g = I. Statistical error bars are smaller than 
symbol size. Dot-dashed lines indicate the HP strong- and 
weak-coupling results. 



in the phonon cloud of an EX-P with center-of-mass mo- 
mentum p. On the other hand, for a given eigenstate 
\v;p) of UJ), the probability W(P,K) for a E-H pair to 
have center-of-mass momentum P and relative momen- 
tum K is obtained by tracing out the phonons from the 
density matrix p — \ v; p)(v; p|: 
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In order to validate this novel technique, we now in- 
vestigate a minimal 3D simple-cubic (tight-binding) two- 
band model, £ c ,„(k) = E C:V ± (E c , v /6) Y, a =x, y ,z( 1 ~ 
cosk a ) , where E c = E g gives the direct gap at k = 0, 
E v = 0, and E c and E v are the bandwidths of the con- 
duction and valence bands, respectively. Furthermore, 
we take the phonon frequency w q = 51, particle-phonon 
couplings <7e,/t(q) = g, and the interband Coulomb at- 
traction C/(p, k, k') = U as momentum independent. 
A distinctive feature of this model is that E and H 
only form a bound state if U > [/», where U" 1 — 
N^ 1 X^keBZ ( £ c(k) - £„(k) - EgA 1 . Hence, depending 
on the Coulomb attraction U and the dimensionless E/H- 
phonon coupling A = 2g 2 /(flE c ), the EX-P binding en- 
ergy Eex-p = E g , s , — 2Ehp (defined as the difference 
between the EX-P ground-state energy E g , s , and twice 
the ground-state energy of a single electron/hole Holstein 
polaron (HP) £hp at p = [l5|]), is either £ex-p = 
(unbound state) or £ex-p < (bound state). 
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FIG. 4: Internal structure of an EX-P in the mass-symmetric 
model with U — 1.5E g and A = 0.195 < A» (1st column), 
A = 0.223 w A» (2nd column), A = 0.29 > A* (3rd column). 
First row shows the integrated Z-factors, second (third) row 
displays the reduced distributions of W^PjK) (O w(K z ) = 
J d s PdK x dK y W (P,K) (w(P z ) = J d 3 KdP x dPyW {P ,K)) . 



Figure [3] presents the data for the mass-symmetric 
model with E c = E v = 3E g and SI = 0.25E g . Here 



the critical Coulomb attraction is [/* ~ 1.98E g , i.e. for 
U/E g = 2.5 the H and E are already bound at A = 0. 
As A is increased, the EX-P binding energy and effective 
mass smoothly increase and finally show the standard 
weak- to strong-coupling crossover [l6| . By contrast, for 
U < £/*, a critical coupling A*(£7) is required to create the 
EX-P bound state. While here the single HP exhibits a 
rather smooth crossover from a weakly to strongly mass- 
renormalized QP at about Awl (see inset in Fig. [3}, 
the transition of the E-H pair from unbound to bound 
state is accompanied by a much more rapid change of 
EX-P properties (Fig. [3|) . The sharper crossover to the 
small radius EX-P regime can be understood from the 
increasing importance of the QP-phonon coupling when 
the EX-QP bound state is established [17| ■ In accordance 



with Ref. we find A*(C7 = 0) » 0.5. 



The internal structure of the p = EX-P state signifi- 
cantly changes at the unbound-to-bound-state transition 
(Fig. 0|). For A < A* the E-H pair is only weakly dressed 
by phonons leading to P ss 0. The distribution of the 
relative momentum K is narrow, pointing to a large real 
space separation of the E-H pair. For A > A*, E and 
H are confined and the K distribution broadens. At the 
same time the P distribution develops a shoulder due to 
fluctuations of the momentum of the phonon cloud, and 
the amplitude of the peak at P = decreases, reflecting 
the suppression of the zero-phonon weight. For A > A* 
the phonon distribution Z^ M ^ is almost Gaussian 14 1. 
In the A w A* region, the Z( M )-factor distribution shows 
a kind of "bimodality" which, however, is only observed 
provided that SI <C E CtV . This feature rapidly disappears 
for larger phonon frequencies [5j . 

For finite SI, the phonon exchange leads to a retarded 
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FIG. 5: (Color online) Exact DMC EX-P energy and effective 
mass (symbols with dotted lines, statistical errors are smaller 
than the symbol size) compared to results for the EAHM 
(solid lines). Data obtained for U — E g and E c = E v = 3E g . 
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becomes almost immobile and the model is equivalent to 
that of a phonon-assisted localization on an attractive 
impurity [13]. The excellent agreement of our data with 
those obtained for the impurity problem by the Cheby- 
shev space method and DMC in direct space [13] pro- 
vides one more successful test for the momentum space 
DMC technique (see Fig. [6]) . Surprisingly, the phase di- 
agram of asymmetric EX-P binding well coincides with 
that of trapping by impurity already at r/ = 0.1. 

In conclusion, we have developed an exact DMC al- 
gorithm for the interacting electron-hole-phonon system 
which fully takes into account an internal structure of 
the exciton-polaron. The technique is applicable to any 
band structure in D—l, 2, and 3, an arbitrary (attrac- 
tive) Coulomb interaction, and momentum-dependent 
fermion-phonon coupling. We constructed the phase di- 
agram for a phonon-assisted electron-hole bound state 
formation, and discussed the change of the internal struc- 
ture of an exciton-polaron in the transition regime. Com- 
parison our exact results for the exciton-polaron problem 
with findings for an effective attractive Hubbard model 
shows that for realistic values of phonon frequencies the 
retardation effects can not be neglected. 

We appreciate helpful discussions with A. Alvermann, 
F.X. Bronold, and M. Hohenadler. E.B. acknowledges 
financial support by DFG through SFB 652. A.S.M. is 
supported by RFBR 07-02-00067a. 



FIG. 6: (Color online) Phase diagram of the mass-asymmetric 
EX-P model with E c = 3E g , E v = r)E c , and fi = E g . Solid 
line (solid circles) indicate transition line of the static im- 
purity problem with r\ = [l8j] (obtained by exact numeric 
techniques pH |20|]). Dashed lines are guides to the eye. 



interaction between E and H. In the antiadiabatic limit, 
E c ,v, the retardation effects become negligible, and 
our model is equivalent to an effective attractive Hubbard 
model (EAHM) with U eS = U + 2g 2 /Q and hopping in- 
verse proportional to the single HP mass t e ff — l/ m HP 
(see, e.g., Ref. [HI). Comparing the EAHM data with 
results obtained for the full model where the retardation 
effects were included provides a good check for our DMC 
algorithm. Indeed we found good agreement for large 
phonon frequencies (see Fig. [5] for f2 = 12E g = 4E C ). 
However, the domain of validity of the instantaneous ap- 
proximation is rather limited. As seen from Fig. the 
effective mass, e.g., considerably deviates from the exact 
result for £1 ~ E c already. 

To make contact with the situations for realistic semi- 
conductors we construct the ground-state phase diagram 
for a set of mass asymmetries 77 = E v /E c < 1 (Fig. [6J. 
The data clearly show that the larger the bare mass of 
the hole m° ~ the smaller U and A are required 

to bind the excitonic QP. For 77 — > (mi — > 00) the hole 
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